Microscopic nonequilibrium theory of double-barrier Josephson junctions 
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We study nonequilibrium charge transport in a double-barrier Josephson junction, including non- 
stationary phenomena, using the time-dependent quasiclassical Keldysh Green's function formalism. 
We supplement the kinetic equations by appropriate time-dependent boundary conditions and solve 
the time-dependent problem in a number of regimes. From the solutions, current-voltage charac- 
teristics are derived. It is understood why the quasiparticle current can show excess current as well 
as deficit current and how the subgap conductance behaves as function of junction parameters. A 
time-dependent nonequilibrium contribution to the distribution function is found to cause a non- 
zero averaged supercurrent even in the presence of an applied voltage. Energy relaxation due to 
inelastic scattering in the interlayer has a prominent role in determining the transport properties 
of double-barrier junctions. Actual inelastic scattering parameters are derived from experiments. 
It is shown as an application of the microscopic model, how the nature of the intrinsic shunt in 
double-barrier junctions can be explained in terms of energy relaxation and the opening of Andreev 
channels. 



PACS numbers: 74.25.Fy, 74.40.+k, 74.45. +c, 74.50. +r 



I. INTRODUCTION 



The Josephson effect is a hallmark of superconductiv- 
ity. It is also the basis of a wide range of applications such 
as metrology, sensing and classical and quantum logic 
circuits. Moreover, the detailed study of the Josephson 
effect in mesoscopic devices provides deep insight in the 
mechanism of the formation and transport of supercon- 
ducting correlations in weak links and at interfaces. It 
is known from the microscopic theory of superconductiv- 
ity that the supercurrent across weak links is carried by 
Andreev bound states (ABS). The supercurrent depends 
both on the ABS energy levels as well as on their pop- 
ulation, i.e. on the quasiparticle distribution function 
over energy. This provides a possibility for external con- 
trol of current. It was realized long time ago that devia- 
tions from equilibrium may strongly modify the transport 
properties of weak links and Josephson tunnel junctions. 
A review of early work on various aspects of nonequilib- 
rium superconductivity is given by the articles in Refs. 
0] and U and by Kopnin^ Two main topics in this field 
are the effects arising from a charge imbalance 4 * 5 " and ef- 
fects from a stimulation of superconductivity by external 
fields^ 

The recent progress in the fabrication of supercon- 
ducting structures of sub-micrometer size has stimu- 
lated a renewed interest in nonequilibrium effects in 
Josephson junctions. The effect of supercurrent con- 
trol by current injection from additional terminals 



was first studied theoreticall y 7 ! 8 ! 9 and demonstrated 
experimentally 10,11 for a diffusive Superconductor- 
Normal metal-Superconductor (SNS) junction. In this 
case, even a sign reversal of the critical current is 
possibleAii Additionally, control of supercurrent by cur- 
rent injection was studied in structures with ballistic 
transport, 12 ! 13 ! 14 ! 15 ! 16 ! 17 

Deviations from equilibrium are also reflected in the 
quasiparticle current. The dissipative current component 
in SNS junctions arises from Multiple Andreev Reflec- 
tions (MAR) of quasiparticles. A quasiparticle gains eV 
in energy each time it traverses the interlayer, resulting 
in a strong nonequilibrium distribution function at sub- 
gap energies, as observed by Pierre et ali& In point con- 
tacts, the microcscopic description of the current in terms 
of MAR was derived by Averin and Bardasji 9 - based on 
a scattering matrix approach, while Cuevas et alm^ de- 
scribed MAR in quantum point contacts by means of a 
tunnel Hamiltonian approach. For the case of an SNS 
junction with a long interlayer as compared to the coher- 
ence length (incoherent regime), the current due to MAR 
was calculated by Bezuglyi et alm^ 

When tunnel barriers (I) are introduced at the SN in- 
terfaces, the quasiparticles in short junctions can undergo 
transmission resonances, resulting in a dephasing of the 
electrons and holes. However, it was shown in Ref. |52 
that for a broad transmission resonance in a SINIS junc- 
tion, the resonance energy width being larger than the 
ABS, coherent transport occurs, and a microscopic model 
was given in terms of MAR, integrated by a universal 
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distribution of transparency eigenvalues^ 3 - It is shown 
by Naveh et al. 24 that the same distribution function 
also describes electrical transport in high critical current 
density junctions. In order to model the nonstationary 
and noncquilibrium transport through SINIS structures 
in the general case, a full Keldysh Green's function ap- 
proach is required. The derivation of this microscopic 
model as well as its solutions is the scope of this article. 

SINIS junctions are promising basic elements for ap- 
plications in classical computing and metrology, because 
they are intrinsically shunted. Both the Rapid Single 
Flux Quantum logio^ 5 - and the digital voltage standards 
are electronic applications for which the use of these 
structures seems promising. 26 However, important fea- 
tures of the IV characterstics are not sufficiently under- 
stood yet, such as the magnitude of the subgap conduc- 
tance and the nature of the intrinsic shunt of these junc- 
tions. The experimental observation that the hysteresis 
in the IV curves depends nonmonotonically on critical 
current density has not yet been explainediS 6 - 2 ^ Thus, 
understanding the transport properties on a microscopic 
level will be valuable for electronic applications. Here, 
these transport phenomena will be clarified in terms of 
the opening of Andreev channels and inelastic scattering 
in the interlayer. 

Earlier work concentrated on modeling the IV char- 
acteristics of double-barrier junctions in specific limiting 
cases. When one of the electrodes is replaced by a nor- 
mal metal, the time-dependencies simplify considerably, 
since formally one can then put voltage to zero in the 
superconductor. The IV characteristics of SININ junc- 
tions were studied by means of the quasiclassical Green's 
functions technique by Zaitsev, 28 Volkov et al.^ and Za- 
itsev et alm& Lempitski: 31 studied noncquilibrium effects 
on the nonstationary properties of long SNS junctions 
in the absence of interface barriers and Kadin~ used a 
time-dependent Ginzburg-Landau approach, valid only 
in a narrow temperature range. Another limiting case 
is the double-barrier structure with a long interlayer as 
compared to the coherence length. The derivation of 
time-dependent transport properties in this case simpli- 
fies since a decoupling of the electrodes is possible, as e.g. 
studied by Volkov and Klapwijk 33 and Bezuglyi et alm^ 

In this article, a microscopic quasiclassical theory will 
be given for double-barrier Joscphson junction with two 
superconducting electrodes and a short interlayer. The 
interlayer will be assumed to be a diffusive normal metal, 
but it will be indicated how the model can be extended in 
a straightforward way to incorporate a superconducting 
gap in the interlayer. The Keldysh formalism is intro- 
duced in section^ The spectral supercurrent density is 
obtained and appropriate time-dependent boundary con- 
ditions are derived to supplement the kinetic equations 
for the energy distribution functions in the interlayer. 
The technical scheme for solving the time-dependent 
Keldysh-Usadel equation may have applications beyond 
the present paper. Solutions are presented in section 
IIIII for the adiabatic limit of eV <C A#. As an intrigu- 



ing nonequilibrium effect in a double-barrier Josephson 
junction, we show that even at finite voltage bias, there 
can be a nonzero averaged supercurrent. Energy relax- 
ation due to inelastic scattering is a phenomenon that 
strongly modifies the energy distribution function. It will 
be shown in section II VI that inelastic scattering is im- 
portant in a double-barrier Josephson junction and how 
this effect can be incorporated in the microscopic model. 
As an application of the microscopic model, the nature 
of the intrinsic shunt of double-barrier junctions will be 
discussed in section The observed nonmonotonic hys- 
teresis vs. critical current density dependence, as well as 
the actual values, are explained. 



II. KELDYSH FORMULATION 

The Matsubara Green's function technique can be ap- 
plied to a many-body system in equilibrium, from which 
the energy-dependent properties of the system can be 
derived. In addition to obtaining spectral quantities, 
we need to know how the states are populated under 
nonequilibrium conditions. For this purpose Keldysh 34 
proposed a set of propagators along a contour in the 
complex-time plane that allows to describe the real-time 
evolution of a system outside equilibrium and at a fi- 
nite temperature. The review of Rammer and Smith 35 
describes the use of the Keldysh technique in the trans- 
port theory of metals. The Keldysh method is introduced 
specifically for nonequilibrium superconductivity in Refs. 

The quasiclassical approximation is used, in the sense 
that rapid oscillations of the wavefunctions on the scale 
of the Fermi- wavelengt are averaged out. Furthermore, in 
this paper it is assumed that the transport through the 
interlayer is diffusive, the thickness being much larger 
than the elastic scattering length, so that the Usadel 
equation can be used. 



A. Time-dependent Usadel equation 

A compact notation of the equations for the quasiclas- 
sical Green's functions becomes possible by introducing 
the Green's function in Keldysh x Nambu space 

The quasiclassical Green's function G is a function of two 
times, t and t' , and the time-dependent Usadel equation 
in the absence of a vector potential reads 3 ^ 

/ — — \ BC BC ~ — 

-MVfGoVG) + r 3 h— + — hr 3 -iA(t)G 

+GiA(t') = -i(£ incl oG-GoE incl ) , (2) 
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where 



B. Retarded and advanced propagators 
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T> is the diffusion constant, Si ne i the self energy with 
retarded, advanced and Keldysh components, * denotes 
the complex conjugate and o denotes a convolution over 

the internal time coordinates, e.g. Si n ci(i, t') o G — 

J dtiSinei (t, ti) G (t%, t 1 ). The function G is normalized 

as G o G — 1 . The expression for the current in the 
Keldysh formalism is 
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The Green's functions can be transformed to energy- 
frequency space (E, us) by Fourier transforming the func- 
tions G(t-t',(t + t') /2), 

G(E,us) = J G (t - t', i±£^ e -iE(t-t')/n 

e^ t+t 'V 2h d(t-t')d(t + t')/2. (5) 

This transformation is analogous to the Wigner repre- 
sentation of the full double- coordinate Green's function. 
Spectral quantities that only depend on energy and not 
on frequency after Fourier transforming, such as the 
equilibrium Green's functions in the electrodes, only de- 
pend on the time difference before Fourier transforming. 
Each term in Eq. J2J can be transformed to (E, w)-space. 
Hence, the Usadel equation can be rewritten in (E,us)~ 
space as 



-DfiV (G o VGJ + iE ^r 3 , oG 
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Equation © consists of an Usadel equation for the 
retarded Green's function, the advanced Green's func- 
tion, and an equation containing the Keldysh Green's 
function. The Usadel equation for the retarded Green's 
function G R in the interlayer (taking the limit of A = 
and zero inelastic scattering, Ej ne ] = 0) in Fourier com- 
ponents reads 

-Pw(G fl oVG B ) +inuso/2{r 3 ,G R (E)} 
+iE[f 3 ,G R (E)]=0, (8) 

where the index n denotes the n-th harmonic. The self- 
energy terms in Eq. JJJJ can effectively be represented by 
a characteristic inelastic scattering time Tj„, as was de- 
rived by Larkin and Ovchinnikov. Taking the inelastic 
scattering to be time- independent is of course a rough ap- 
proximation, but suffices for the purposes of this paper. 
The self-energy terms have been neglected in Eq. (jHJ, 
which is justified as long as H/Ti n <C fcsT. Note, that 
Eq. JHJ can be reduced to the time-independent case for 
n = 0. Time-dependence occurs if the boundary con- 
ditions, which will be detailed below, provide nonzero 
limiting values for G n . 

The most general decomposition of the retarded and 
advanced Green's functions is a linear combination of the 
three Pauli matrices^ Whenever the phase is constant, 
it can be chosen such that it suffices to define 



qR{A) =G n(A) h+F R(A)~ u 

G A = -fsG^f, = -G R *~ 



T3 



F R *f v 



(9) 



Assuming that the thickness of the interlayer d is much 
smaller than the coherence length in the interlayer, £ = 
y/T>/ ' E where E is a characteristic energy at which the 
system is probed, we can take the retarded Green's func- 
tion much larger than its gradient. The double-barrier 
structure under consideration is depicted in Fig. ^ Inte- 
grating both sides of Eq. (JSJ over the interlayer thickness 
and barriers gives 



where A (f) = is taken for simplicity. [t 3 ,G] is 

the commutator of r 3 and G, and {r 3 ,G} is the anti- 
commutator. A decomposition of the Green's functions 
in Fourier harmonics can formerly be introduced as 



G(E,u)= G n (E)6{us 



u ) 



(7) 



where G n [E) = G (E,nu>o), as was for example done in 
Ref 20. The generation of higher order Fourier harmon- 
ics is a manifestation of the nonlinearity of the device, 
prevalent e.g. in Eq. ©. 
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+m^-d {f 3 , G R (E)} + lEd [t 3i G r (E)] = 0.(10) 

Zaitse\i2i derived effective boundary conditions for the 
quasiclassical Green's function formalism. These were 
further developed for diffusive scattering in the interlayer 
by Kupriyanov and Lukichev.- 2 Using the Kupriyanov- 
Lukichev boundary conditions for the retarded Green's 
functions. 
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FIG. 1: Schematic representation of the double-barrier SINIS 
structure. Two superconducting electrodes (S) are separated 
by two delta-shaped potential barriers (Ii and I2) and a nor- 
mal metal interlayer (N). The position dependence of the pair 
potential has been indicated by shading. 



where 73 = R B /p£,, p is the resistivity of the interlayer, 
and R B is the interface resistance, we obtain 



+2irk B T cS 



T 3 ,G* (E)\ lB d/(, 

= 0, (12) 



where G R are retarded functions in the left and right 
electrodes respectively. The normalization condition for 
G R in energy-space and decomposed into Fourier har- 
monics can be found from the expressions of Appendix 
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(13) 

Equations (|12() and (|13|) form a complete set of equations 
from which the Fourier components of G R can in prin- 
ciple be determined. This recursive schemes reflects the 
fact, that superconducting correlations can be induced 
over several MAR cycles. Solving the set of equations 
is complicated by the recurrent nature of the equations. 
The Fourier harmonics are coupled to each other and 
have arguments that are shifted in energy. 

From the full set of equations, we can find back the 
quasi-stationary Matsubara case by keeping only the 
n = harmonic of G R and the n = ±1 harmonics of 
F R and neglecting energy-shifts in the arguments. This 
provides a solution for G R and F R that coincides with 
the analytical continuation (oj — » — iE) of the Matsubara 
solution at ip = ir/2. The general Matsubara solutions 
for double-barrier junctions with 751,2 ^ 1 and g?/£ -C 1 
were obtained in Ref. |2(| In the limit of A = the an- 
alytical continuation uj —> —iE of this solution provides 
the Green's functions at T — as function of energy 
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FIG. 2: Normalized density of states in the interlayer at 
T = for several values of the suppression parameter 7 e fl> 
The inset shows the minigap that is present for j e g = 10 2 on 
a smaller scale. 



where F s = A s / (A| - E 2 ) and G s = Ej (E 2 
The asymmetry and effective suppression parameter are 
respectively 
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(15) 



As an illustration to these retarded Green's functions, 
the density of states in the interlayer, N = ReG, is shown 
in Fig. |21 It can be seen that for 7 e g ~S> 1 the den- 
sity of states is determined by a minigap with a value of 
cos((^/2)7rfcsT c /7eff . In the coherent regime^ of y e g -C 1 
the gap in the density of states is given by Acos(y>/2). 
The density of states for intermediate values of the sup- 
pression parameter is characterized by a two-peak struc- 
ture. These findings coincide with the calculations of 
Bezuglyi et ulJ^ in the limiting case of a short interlayer. 



C. Spectral supercurrent 

Supercurrent is carried by states in the weak link and 
their occupation is determined by a distribution function. 
The supercurrent-carrying density of states, or spectral 
supercurrent Im/5 (E), can be determined from G R and 
G A by 
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(16) 



The supercurrent in the regime of 7 c ff <C 1, is found to 
have a spectral density 
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FIG. 3: Normalized spectral supercurrent density as a func- 
tion of energy for various values of the suppression parameter 
7eff. The phase difference between the superconducting elec- 
trodes was fixed at (p = n/2. The inset shows the spectral 
supercurrent in the coherent regime of 7 c ff <C 1. 



for A s cos(c/?/2) < E < A s , while luiI s (E) = for 
E < As cos((p/2) and E > As- This universal ex- 
pression is independent of the interlayer thickness, bar- 
rier height and contact dimensionality as long as the 
number of conduction channels is large4£ The same ex- 
pression was found in the case of ballistic interlayer 
transport^ The spectral supercurrent is nonzero only 
in the range Agcos(i^/2) < E < As, i.e. there is a 
minigap As cos(</?/2) in the spectrum of the Andreev 
bound states, see the inset of Fig- EI O n the other hand, 
all states in the energy range As cos((/?/2) contribute to 
the supercurrent. In long junctions, similar behavior for 
the DOS was found in Ref. |4J and for the current in 
Ref. |a The contact is in the intermediate regime be- 
tween a short ballistic SNS weak link, with bound state 
energy As cos((/?/2) and a tunnel junction with bound 
state energy As- Physically this is caused by the prop- 
erties of the distribution of transparencies, which is a 
combination of open and closed channels (see Ref. l22j) . 

In the incoherent regime of 7 e fi 3> 1, this universality 
breaks down. The minigap in the spectrum of Andreev 
bound states is now given by cos(ip/2)TrkBT c /"f e g. Fig- 
ure |3] shows the spectral supercurrent density for several 
values of the suppression parameter. The sign change 
at E = As has also been observed by Bezuglyi et al~ 
and Heikkila et alm& Going beyond the approximations 
d <C £ and jb *C 1, it was shown by Schapers et aZ41 
for ballistic junctions that low-energy states are gradu- 
ally filled in for larger interlayer thickness and by larger 
barrier transparency. 



D. Kinetic equations 

The energy distribution functions, that determine the 
occupation of spectral funtions, can be determined from 
the kinetic equations. 

From the matrix normalization condition G o G = 1 , 
the upper right component implies that G R o G K + G K o 
G A = 0. Hence, G K can be parametrized as 



G 



K 



G R 0/-/0G' 



(18) 



Furthermore, it was shown by Schmid and Sch6na& and 
Larkin and Ovchinnikov, 40 that / can be chosen to be 
diagonal. We will adopt the notation 



f = f L l + f T T 3 , 



(19) 



where fa and fa are those parts of the distribution func- 
tion that are respectively even and odd in energy. There- 
fore they are named longitudinal and transverse energy 
distribution function respectively. The functions can be 
identified with energy and particle flow42 Physically, a 
deviation of fa from equilibrium is associated with a dif- 
ferent effective temperature and a deviation of fa from 
equilibrium with a chemical potential shift. In equilib- 
rium, fao = and /lo = tanh(£ , /2fcsT). 

Putting Eqs. (|18|) and l|19fl into the Keldysh component 
of the Usadel Eq. J2J an d by making use of the Usadel 
equations for the retarded and advanced Green's func- 
tion, finally the kinetic equations for the Fourier compo- 
nents of fa and fa can be written as 

(D L o V 2 /l)„ + (ImJs o V/ T )„ = (moo + hrr 1 ) 
x [G R o (fa - f 6 n0 ) - (fa - JbM G A ] n , 

(D T o V 2 /t)„ + (Im/s o Vfa) n = -L (niu; + hr^) 
x(G R ofa~faoG A ) n , (20) 

with the generalized transverse and longitudinal diffusion 
coefficients being AD T = Tr(l - f 3 G R o f 3 G A ), AD L = 
Tr(l — G R o G A ). With the parametrization of Eq. © 
this can be further rewritten as Dt — (ReG) 2 + (ReF) 2 
and D L = (ReG) 2 - (ImF) 2 . In obtaining Eqs. (j2U|) . 

use has been made of the rewriting of the S o G — G o E 
term by Larkin and Ovchinnikov^ into a collision inte- 
gral with characteristic inelastic scattering time Ti n . A 
has been assumed to be negligible for simplicity, but a su- 
perconducting gap in the interlayer can be incorporated 
in the model in a straightforward way by keeping the 
terms in the Usadel Eq. J2J that depend on A. In the 
limit of slow time variations, a Fourier transform over the 
time difference provides the known mixed representation 
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of the kinetic equations 38 

DD L V 2 / L +PIm/ s V/ T 

-ReG (V + (f L - /o) = 0, 

VD T V 2 f T + VlmI s Vf L 



-ReG(^ + -l/ T = 0, 



(21) 



which follows directly from Eq. (|2U[1 for the lowest Fourier 
harmonic. The expressions for the supercurrent and dis- 
sipative current components can be derived^ from Eq. 

m 



Is = 



2eR N 



dEf L (E) Imls(E), 



In 



J dED T (E) V/t (E). 



(22) 
(23) 



What remains to be derived, is a proper set of time- 
dependent boundary conditions for the kinetic equations. 



E. Time-dependent boundary conditions 

The Kupriyanov-Lukichev boundary conditions^ for 
the quasiclassical Green's functions can in general be 
written as 

7s£Go^G = GoGi-GioG, (24) 
dx 

where Gi and G denote the Green's functions at the two 
sides of the first interface. From the definition of the 
Green's functions in Keldysh space, Eq. QJ, a boundary 
condition can be written for each matrix element. In 
Appendix [XI this set of boundary conditions is rewritten 
into 



7b£ 



1 - G R G A 



d -jL+(r,-G R f,G A )4zh 



dx V I dx 

G R (Gf - Gf ) - (Gf - G A ) G A ] (f L1 - f L ] 

G R (Gf fa - r 3 Gf) - (Gf fa - f 3 Gt) G A 



x {hi - /t) , 



(25) 



where all the products have to be regarded as time- 
convolutions and /li/ti are the distribution functions 
in the respective reservoir. 

The Green's functions become time dependent by ap- 
plying a voltage over the interface. In the absence of volt- 
age, the Green's functions in the electrodes only depend 
on time difference since equilibrium is assumed. The po- 
tential can be introduced in each electrode by a Gauge 
transformation 41 of the Green's function in the electrodes 



Gf (A) (t, = S (t) G R[A) (t - t') £t (t>) 



(26) 



where S (t) and Sft (t r ) are given by 



S(t) 



St ( t ') = 



ieVt/h 

( 

-ieVt' /h 







-ieVt/h J 1 





JeVt'/H 



(27) 



Volkov and Klapwijk 33 performed a Gauge transforma- 
tion of the interlayer Green's functions, which works only 
in the the limit d 3> £ because of the small coupling be- 
tween the electrodes in this case which allows to neglect 
the interference terms in the interlayer leading to a local 
time-dependence. 

By performing the Gauge transformation for Gf and 
G A and by taking the trace from Eq. I|25l) . one obtains 
the first boundary condition in time representation. The 
second equation is obtained by taking the trace after mul- 
tiplying left- and right-hand side of Eq. (|25|) by f a . This 
results in 



DrlBi-rfr = { ReGiReGisin 
dx ' 



eV 



(t-f) 



eV 

+ImiqReFsin — 
n 

— fx jReGiReGcos 
+ReFiReF cos 



(t + tf) 



(/lo - /l) 



eV 



(t-f) 



T (t + f) 



(28) 



DlIb^Il = { ReGiReG cos 



TmiqlmF cos 



eV 



-fx |ReGiReGi sin 
~eV 

KeFilmF sin 



(t + O 

eV 

IT 

(t + 



(/lo - fh) 



(t-f) 



(29) 



where all products are time convolutions and use has 
been made of the fact that Jti = 0, since the electrodes 
are assumed to be in internal equilibrium. The energy 
distribution functions are not only coupled through the 
kinetic equations l|2(JI) . but through the boundary condi- 
tions as well. 

At the second interface a similar set of boundary con- 
ditions can be derived, which can be obtained from Eqs. 

and J23) by replacing Gi and Fi by G 2 and F 2 re- 
spectively, and by multiplying the right-hand side of Eqs. 
© and (US) by -1. 

Note, that D L = (ReG) 2 - (ImF) 2 = for energies 
smaller than the minigap in the interlayer. Hence, for 
energies at which Dl = 0, boundary condition (|29|l is 
replaced by fi, = /lo- This physically means, that the 
system does not conduct heat inside the gap and that the 
distribution in the gap is controlled by coupling to some 
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external heat bath, e.g. through the substrate, and not 
through the superconducting leads. 

Each term in the boundary condition contains time 
convolutions. With the aid of the expansion of the 
Green's functions in Fourier harmonics and the expres- 
sions of Appendix iBl for the time convolutions of double 
and triple products, the convolutions can be worked out 
for each term. The left-hand side of Eq. (|28[l is for ex- 
ample 

oo 

d f 

Dt ° 1b ^^ t = Y1 I D T ,n(E + n'Lj () /2) lB Z 



dx 



(30) 



The sine and cosine dependencies in the boundary con- 
ditions cause additional voltage shifts as well as coupling 
to higher harmonics, which can be seen for example in 
the term 

f L o [ReG s o ReG o i sin (t - t')] 

= J2j dE fL,n (E + n'LO /2) RcG„' (E - nu /2) 

n.n' 

iE'(t-t')/% e imfu> {t+t') 



xe 
1 
2 



n' — n — 1 
ReG s E + luq 



(31) 



In principle, the set of kinetic equations (|2U[I together 
with the boundary conditions Eqs. I|28[) and l|29[l . and 
expressions for the time convolutions, such as Eq. i|30|) 
and Eq. (S3J, now provide a complete set of equations 
to solve the energy distribution functions as function of 
voltage. However, the coupling to higher harmonics and 
energy shifts within the functions themselves make solv- 
ing the equations cumbersome. In principle, a solution 
should crossover to the solutions as found by the MAR 
approach^ in the limit of 7 e g C 1. In the next section 
an adiabatic approximation will be developed in order to 
solve the kinetic equations for eV <C As and a larger 
suppression parameter. 



III. ADIABATIC DYNAMICS 

A. Adiabatic approximation 

In order to simplify the time-dependencies, an adi- 
abatic approximation can be made. When voltage is 
small, the phase oscillates slowly and can even be con- 
sidered quasi-stationary. In this case, we only need to 
keep the time dependence in expressions that contain the 
phase, but can neglect all other time dependencies. Con- 
sequently, the time convolutions become simple products 



and the energy shifts can be neglected. Therefore, this 
approximation is called adiabatic. 

A formal derivation of the parameter regime in which 
the adiabatic approximation can be used, is based on 
the time dependence in Eqs. lj2"5|) and The 
quasiparticle current is determined by the left-hand 
side of Eq. (|28p. namely -Dt7b£^/t /dx. It will be 
shown in this section that the right-hand side of this 
boundary condition in the adiabatic limit is equal to 
the Jt terms in Eq. (|28|) . Hence, deviations from 
the adiabatic approximation in the quasiparticle cur- 
rent are only to be expected when the terms propor- 
tional to in Eq. (|28l) are not negligible. The first 
of these terms is fi,ReGiReGi sin[eV (t — t')/h], which 
can be neglected for eV <C As- The second term is 
fi,lmFiReF sin[eV(t + t')/H), which is nonzero only due 
to the loop- like construction with Eq. (|29|l . in which 
the terms ImFi and ReF are shifted eV/2 in energy ev- 
ery cycle, making their overlap nonzero after approxi- 
mately As/eV cycles. For large suppression parameters, 
ReF ~ 7 cff 1 . Hence, smallness of this term can now be 
formulated as (l/"/ c g) As ^ eV <C 1. Combining the con- 
ditions, the conclusion is reached that the adiabatic ap- 
proximation is valid when eV <C As and j c ff 1. 

In this case, the phases xi,2 can be introduced by the 
parameterization 



W,2 



G? 2 F&e***.* 



P« p-*Xl,2 



-G? 



(32) 



No additional Gauge transformations have to be per- 
formed in the boundary conditions. Hence, we can use 
the parametrization of Eq. (|32|l directly in the bound- 
ary condition (|25fl . The first boundary condition is then 
found by taking the trace of Eq. (|25|l . The second is 
found by taking the trace after multiplying with T3. With 
Eq. (l.')2[l and after some rewriting this gives 



lBi^D T — f T (0, d) = ±M T1>2 [f T (0, d) t /to] , 

lBi,^D L ^f L (0, d) = ±M Lh2 [f L (0, d) - f L0 ] , (33) 
ax 

where 

1 E + eV/2 1 E-eV/2 . . 
fw.To = 2 tanh ± - tanh , (34) 

arc the distribution functions in the leads and 
M T1 , 2 = ReGs^ReG + ReFsReFcos( X i - X2), 
Mli,2 = ReGs^.^ReG — ImFsImF cos(xi — X2), where 
Xi — X2 = v/2 = eVt in the case of symmetric barriers. 
Here, G and F are given by Eq. I|14(l . and 



E±eV/2 _ As 

(F± e y/2) 2 -A| V^s - E2 



(35) 



As can be seen from here, using also section III Dl the 
kinetic equations in the quasi-stationary limit coincide 
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with the known adiabatic equations in the mixed rep- 
resentation as derived by Larkin and Ovchinnikov^ In 
the absence of inelastic scattering, the equations further 
simplify to 



0.20 



D T {E) 
Dl (E) 



d 2 f T 
dx 2 
d 2 f L 



(E,x)+ImIs(E)^(E,x) =0, 
ax 



(j)2 {E,x)+ImI s (E)-^(E,x) = 0. (36) 

Therefore, Eqs. (|33|1 and l|36|l provide the set of equa- 
tions to describe the time-dependent transport in double- 
barrier junctions in the limit of a large suppression pa- 
rameter and a small voltage. 



B. Dissipative current 

From the expression for the dissipative current compo- 
nent in Eq. (|23() . and by solving Eqs. (|33|) and l|36|) with 
an Ansatz /t,l = ai,2^ 2 + ^i,2^ + ci,2, it is obtained that 



hv (*) 



1 



cRn 



dE J B i + 7B2 



7-Bl 

Mti 



1B1 

M T2 



/TO- 



(37) 



This solution is obtained in the limit of no inelastic scat- 
tering. The effects of energy relaxation in the interlayer 
will be discussed in section lTVBI It follows from Eq. I|37|) 
that the quasiparticle current has in general a phase- 
dependent contribution through the coefficients Mtx,2- 
The current is therefore time-dependent since the phase 
difference is given by ip — 2eVt/h. The dc component is 
then determined by averaging over time. 

When one of the superconducting electrodes is replaced 
by a normal metal, the expressions for Mt,l simplify, 
since Fn — 0. By putting voltage to zero in the supercon- 
ductor, it can be shown that the expression for the quasi- 
particle current, Eq. 137fl with Mt2 =RcG, coincides 
with the known results of the SININ junction of Volkov 
et alM The additional term m(E) = d -1 J dx/M T (E, x) 
in Ref. is neglected in our case, since the term is small 
as compared to 7_b/Mt. 

As a measure of the subgap conductance of a double- 
barrier Josephson junction, the conductance at eV = As 
is calculated and shown in Fig. 0] as function of temper- 
ature in the limit of 7 e s 3> 1. The inset of Fig. 0] shows 
the conductance at eV = Aj as function of the inverse 
suppression parameter. It can be seen that the conduc- 
tance is enhanced by a decrease in j c g. Physically, this 
corresponds to the opening of Andreev channels due to 
the term ReFsReF in My. For 7 c ff > 10, for which 
the model of this section is a good approximation, the 
conductance is found to be approximately proportional 
to . This proportionality will be used in section 
to predict an intrinsic shunt in high-J c double-barrier 
Josephson junctions. 

High voltage bias. In the regime of a large voltage 
eV ^ Ag, the time dependencies in the electrodes be- 
come decoupled and the current in an SINIS junction can 




0.0 0.1 0.2 0.3 0.4 0.5 

k B T/A s 

FIG. 4: Normalized conductance at eV = As as function 
of temperature (normalized to As) for 7 c ft S> 1. The inset 
shows the normalized conductance at eV = As as function of 
l/7eH at k B T/ A s = 0.5. 
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FIG. 5: Excess and deficit current as function of asymmetry 
for several values of the suppression parameter. 



be seen as the summation of the current in an SININ' 
junction and an NTNIS junction. In this case, the rele- 
vant functions Mt in Eq. (|37(l simplify, and the presence 
of excess and deficit current can be calculated. Figure [S] 
shows the resulting dependence of the excess and deficit 
current on the asymmetry parameter, for several values 
of the suppression parameter. Note, that the decoupling 
into an SININ' and NTNIS junction at eV > A s is valid 
for all values of j c g. The limiting case of a deficit current 
eldefR-N = 4As/3 for the symmetric limit and j e s 3> 1 
coincides with the findings of Zaitsev^i and Volkov et 
al^tL The excess current el ex R^ ~ 1.05 As for 7 e g <§; 1 
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coincides with the results of MAR calculations^ 



C. Nonequilibrium supercurrent at finite voltage 

From Eq. 11' Ml and a solution of the kinetic equations, 
the supercurrent can be determined as a function of volt- 
age. In most tunnel junctions and weak links, the time- 
dependence of the spectral supercurrent is harmonic and 
by averaging over time, the supercurrent becomes zero 
at a finite voltage. However, due to the additional time 
dependence of /l, the product of /l and Im/5 not neces- 
sarily has to be harmonic, and a nonzero time-averaged 
supercurrent can exist at a finite voltage. Physically, the 
time-dependence of /l originates from the fact that due 
to the proximity effect, the heat diffusion coefficient de- 
pends on the phase difference, as is known e.g. from 
Andreev interferometers^ We would like to calculate 



0.0015 



Is (V) = Q f L (E, t)lml s (E, t) dE 



(38) 



where the brackets denote time-averaging. It can be seen 
that this expression explicitly depends on the spectral 
supercurrent, which could for example be suppressed by 
a magnetic field. Therefore, this current contribution is 
true supercurrent in the presence of an applied voltage. 

Since the time-dependent perturbation of Jl is much 
smaller than /o, the kinetic equations H21[) can be rewrit- 
ten in the adiabatic form like Eqs. (j2EJ), now including 
inelastic scattering, 



D 



d 2 f L 
dx 2 



d 2 fr 
' dx 2 



Im/i 



dfr 
dx 



dx 



fo), 



(39) 



where, S is introduced as <5 -1 = N£ 2 /T>Ti n . Under the 
conditions of the adiabatic approximation, we can use 
again boundary conditions l|33[) . and with Ansatz solu- 
tions fr — a\x 2 + b\x + c\ and /j, = a2X 2 + b2X + c 2 , this 
provides us with the solution 



Jl - fo = 7cfflm/s 



M 



1bD t 



fro 



(40) 



where 



M 



2M T1 M T2 + 2-f eS S~ 1 M T1 



(Mti + M T2 + 27cff<^ 1 ) (Mil + M L2 + 2 7eff ^- 1 ) ' 

(41) 

From Eqs. (|40|) and 141|) it can be seen that in the limit 
of strong inelastic scattering, Ti„ — > 0, Jl — fo will be 
proportional to Ti n and therefore equilibrium is restored. 

The nonequilibrium correction to Jl is obtained within 
the adiabatic approximation. Therefore, we assume that 
time dependence only comes into the final expressions 
via the phase factor in the spectral supercurrent density 
Im/5. In the case of symmetric barriers, Im 2 /g averaged 
over time is equal to the average of sine-squared, which 
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FIG. 6: Nonequilibrium averaged supercurrent at a finite 
voltage bias (eV — O.lAs) at feT = 0.2As and a small in- 
elastic scattering rate 7 = 10 _3 As as function of the suppres- 
sion parameter. The inset shows the supercurrent as function 
of bias voltage at fcgT = 0.2Ag for a fixed suppression pa- 
rameter 7 e fl. 



is just a factor 1/2. The supercurrent contribution can 
therefore be written as 



1 



2eR 



-7cff 



N 



00 

J D^ 1 (Im/5) 2 f T0 MdE. 



(42) 



In order to perform the integral some smearing of the 
Im/| divergency has to be assumed. Physical reasons 
for this are always present, like a small amount of in- 
elastic scattering. An inelastic scattering term 7 can be 
taken into account in the retarded part of the Usadel 
equations^ but in the limit of little inelastic scattering, 
7 *C ksT c s, the scattering term can be presented in the 
solutions by transforming the energy E to E + 17. In 
section realistic values of the inelastic scattering pa- 
rameter and the suppression parameter will be derived. 
From these values it can be concluded that 7 e ff<5 -1 <C 1 
in most of the practical cases. In this limit, the only 
contribution to the supercurrent comes from E > As, 
the spectral supercurrent being zero below the minigap 
and M being zero between the minigap and A5. The 
resulting supercurrent in this limit is shown in Fig. as 
function of the suppression parameter and in the inset as 
function of voltage. 

The physics of this effect is similar to the physics 
of the nonequilibrium supercurrent first considered by 
Lempitskii 31 for a long diffusive SNS junction without 
potential barriers at the NS interfaces. The mechanism 
is the conversion of quasiparticle current into supercur- 
rent inside the junction, formally described by the cou- 
pling term Imlsd/x/dx in Eq. (|39|l . An alternative ex- 
planation for this mechanism has been given in terms of 
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thermoelectricity in Ref. 0. However, quantitative dif- 
ferences occur between the cases of long and short junc- 
tions. In a long SNS junction with d S> £ the strongest 
deviations of the distribution function Jl from equilib- 
rium occur at the sub-gap energy range, at energies of 
the order of the Thouless energy KD/d 2 . In the case of 
MAR, an additional non-equilibrium correction to /l ap- 
pears, as shown by Pierre et al.^ which is beyond the 
present adiabatic approach in which voltage and interface 
transparencies are small so that MAR is suppressed. At 
sub-gap energies the excitation of the symmetric mode 
described by /l, generated by the quasiparticle current, 
cannot diffuse out of the junction since the corresponding 
diffusion constant Dl vanishes in the S-electrodes. On 
the other hand, in SINIS junctions, the symmetric mode 
at low bias is excited only at E > As since the quasipar- 
ticle current vanishes in the energy range E < As due to 
the presence of tunnel barriers. Hence, deviations of fi 
from equilibrium occur only at E > A5. Since the dif- 
fusion constant in the S-electrodes Dl < 1 at E > As, 
the excitations at E > As are partially trapped at these 
energies, though the magnitude of the effect in SINIS 
junctions is smaller than in SNS junctions. With a de- 
creasing barrier heigth in SINIS junctions, an Andreev 
contribution to the quasiparticle current appears and de- 
viations from equilibrium should also occur at E < As- 
The study of this crossover, however, is beyond the scope 
of this paper. 

The noncquilibrium supercurrent at a certain voltage 
is approximately two order of magnitude smaller than the 
dc supercurrent at zero voltage, but is comparable with 
the quasiparticle current at the same voltage. The latter 
conclusion is based on the assumption that inelastic scat- 
tering is negligible, but we will see in the next section that 
energy relaxation in the interlayer increases the subgap 
quasiparticle current. Therefore we can conclude that the 
nonequilibrium supercurrent can be present in double- 
barrier junctions, but that under realistic junction pa- 
rameter values, the contribution to the total current is 
minor. Note, that the nonequilibrium supercurrent in 
the considered regime is proportional to the square of the 
sine of the phase difference over the junction, Is ~ sin 2 </>. 
This effect can give rise to the occurence of half-integer 
Shapiro steps in the IV characteristics when the junction 
is irradiated by microwaves 1 50 ' 5 1 i 



IV. INELASTIC SCATTERING 

In many mesoscopic systems and weak links, the time 
of flight of a quasiparticle through a normal metal or 
superconducting layer is much shorter than the char- 
acteristic inelastic scattering time in the specific mate- 
rial. Hence, inelastic scattering in mesoscopic systems 
and weak links is usually neglected. However, in double- 
barrier junctions, time of flight can be large. Because 
of the normal reflections at the interfaces, a quasiparticle 
on average traverses the interlayer many times. The time 



that a quasiparticle effectively spends in the interlayer is 
proportional to D , where D is the transparency of each 
barrier. For a transparency of the order of 10 -6 , the time 
of flight in the interlayer is for example of the order of 
t = d/DvF = 0.5 ns, for a thickness of about 10 nm and 
a typical Fermi- velocity of 1.5 x 10 6 m/sS 

In most double-barrier Josephson junctions Al is used 
as an interlayer material and therefore the inelastic scat- 
tering time in Al should be considered. Kaplan et alJ& 
estimated an inelastic scattering time in bulk Al of 400 
ns which is much larger than 0.5 ns. However, magne- 
toresistance and microwave measurements in thin films of 
Aj54j55 gho-^d that the inelastic scattering time in thin 
Al films is orders of magnitudes smaller than in bulk, 
namely of the order of 0.1 to 1.0 ns in films of a few 
to 10 nm thickness. Therefore, in the modeling of time- 
dependent transport properties of double-barrier junc- 
tions, inelastic scattering, or energy relaxation, has to be 
taken into account. The inelastic scattering comprises 
both electron-phonon and electron-electron scattering. 



A. Derivation of a microscopic model 

In this section, a microscopic model will be derived 
for the quasiparticle current as function of voltage in 
double-barrier Josephson junctions with low-transparent 
barriers. It will be shown that the results coincide with 
the phenomenological model by Heslinga and Klapwijk, 56 
who derived their model by matching the population and 
extraction rates of the quasiparticles in the interlayer. 

In this section, the assumption will be made that 
7cff 3> 1- In this approximation, the proximity effect 
can be neglected, i.e. — 0. Furthermore, the 

spectral supercurrent lm!s(E,t) = 0, ReG (E) — 1, and 
Dl = Dt = 1- In this case, none of the quantities ex- 
plicitly depends on time. Then, the kinetic Eqs. (|20[) can 
be simplified to 

V ^n-^h (E, x) - [f L (E, x) - f (E)] = 0, 
'DT in -^f T (E,x)-f T (E,x) =0, (43) 

where /o (E) = tanh(i?/2T) and T> is the diffusion con- 
stant in the interlayer. Use is made of the fact that 
/to (E) = in equilibrium. The kinetic equations are de- 
coupled in this case, but /t and Jl are coupled through 
the boundary conditions. 

The boundary conditions can either be obtained by 
simplifying the relevant terms of the expressions that 
contain all harmonics, such as Eqs. (|3U|) and (|3 II) . or by 
starting from the time-dependent boundary conditions, 
Eqs. (|28[1 and i|29|) . In the latter case, the transforma- 
tion to energy space is straightforward. The right-hand 
sides of Eqs. (|28[1 and (|29|) only contain terms that de- 
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FIG. 7: Semiconductor-diagram representation of tunnel and 
scattering rates in a double-barrier junction at bias-voltage 
eV = As. The energy conserving processes, (i) and (ii), are 
in the case of inelastic scattering complemented by process 
{Hi). 



pend on time difference since F R ^ = 0, e.g. 

eV(t-ti) 



KeG\(t — ti) o isin ■ 



o ReG(ii - t') 



of L (h -t')= J dEe-^-^ReG (E) f L (E) G_, 

(44) 

where G_ = Red (E + eV/2) + Red (E - eV/2). The 
left-hand side of Eq. (|28[) becomes 

o 

DTlB^f T (t,t') 

= J dEe- lE ^ 1 ') D T (E) f T (E). (45) 

Hence, together with RcG = 1, finally the boundary con- 
ditions read 

7s^/t (E, ±d/2) = t/t (E, ±d/2) N + 
-f L (E,±d/2)N_+R_, 

7b^/l (E, ±d/2) = t/l (E, ±d/2) N + 

-f T (E,±d/2)N-±R-, (46) 

where N± = ReGi (E + eV/2) ± ReGi (E - eV/2) in 
the superconductors and R± — ReGi (E + eV/2) x 
f (E + eV/2)±ReG 1 (E - eV/2) f (E - eV/2). The ki- 
netic equations provide that c\ = 2a\Dri n = 2a\8 for the 
Ansatz fx — a\x 2 + b\x + c\ and /j, = a-ix 2 + 62^ + c%. 
Using boundary conditions l|46(l and neglecting terms 
proportional to d 2 , the solution can be simply found. 
The quasiparticle current is given by Eq. (|23|) . where 
dfr/dx — bi, and b± is given by 



1 R-N+ - R+N- + (R_ - f N_) /rr„ 



7b £ 



N 4 



i/rr„ 



(47) 



where T" 1 = j B d£/Vh = e 2 N(0)R B d/h is the tun- 
neling injection rate into the normal metal interlayer, 



iV(0) the unnormalized density of states in the inter- 
layer and Rb the specific barrier resistance. Using the 
fact that density of states functions are symmetric in en- 
ergy and /o is asymmetric in energy, (i?_ — /o-/V_) can 
be simplified to 2ReG 1 (E + eV/2)[f (E + eV/2)-f (E)] 
and R-N+ - R+N- = 2ReGi(£ + eU/2)ReGi(£; - 
eV/2)[f (E + eV/2) - f ((E - eV/2))}. With a N /2 lB = 
e 2 N(Q)V/ lB = R~\ the expression for the quasiparticle 
current finally becomes 



eR 



N 



00 

J dEReGi {E + eV/2) 



ReGi (E -z£)F_+ [f (E + zf)- f (E)] /T n , 



i/rr„ 



(48) 



where F_ = f (E + <%-) - f (E - ^-) and G+ = 
ReGi (E + eV/2) + ReGi (E - eV/2). This expression 
is equivalent to the findings of Heslinga and Klapwijk, 56 
in the limit of RcG = 1 , who derived a model by equat- 
ing the population and extraction rates in the interlayer. 
Zaitsev's results for the SININ junction in the limit of 
no energy relaxation^ coincide with our findings as well. 
The equivalence of a mesoscopic or phenomenological ap- 
proach and the more rigorous Green's functions treat- 
ment is shown by ArgamanSi to hold for the equations for 
current. Here, we have proven that the final expression, 
Eq. (|48|l also follows from the Green's function approach, 
using the appropriate boundary conditions. 



B. Influence of inelastic scattering on transport 
properties 

Examples of possible tunneling processes are indicated 
in Fig. [7J In one of the processes a quasiparticle is in- 
elastically scattered in the interlayer. Equation l|48|l coin- 
cides in the limit of strong inelastic scattering (TTi n = 0) 
with the known result for two SIN tunnel junctions in 
series. In the absence of inelastic scattering, Eq. (|48|) 
reduces to 



2dE ( eV\ f eV\ F_ 

—ReGi (E + -^ Gl (£- T )g-. 

(49) 



Figure |H1 shows both limiting cases as well as IV curves 
for intermediate values of the scattering parameter, taken 
in the present limit of 7 ff 3> 1. It can be seen in the inset 
of Fig. |H1 that inelastic scattering enhances the subgap- 
conductance. This effect will be discussed in section Ivl in 
order to explain the large subgap conductance observed 
in double-barrier junction measurements. 

Equation (|49|l gives a deficit current of elde/RN = 
4A5/3 for eV 3> A5, which coincides with the findings 
of section IIII BI In analogy with the approach of section 
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FIG. 8: IV characteristics at k B T/A s = 0.25 and 7 cff » 1 
on the basis of Eq. 1481 for several values of the inelastic 
scattering parameter IYm. The inset shows the subgap con- 
ductance at eV — A as function of Tn n . 
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FIG. 9: Excess and deficit current as function of the supres- 
sion parameter for several values of the inelastic scattering 
parameter Yn n . 



1111 Bl to calculate excess and deficit currents by summing 
the respective contributions from SININ and NINIS junc- 
tions, the same can be calculated by including inelastic 
scattering as well. Figure [5] shows the resulting crossover 
from excess to deficit current as function of the suppres- 
sion parameter for several values of the inelastic scat- 
tering parameter. For rr, n < 10 _1 only a small deficit 
current is predicted at eV ^> As- However, at moderate 
values of V, i.e. 2A5 < eV < 4Ag, still a considerable 
deficit current is present, as for example can be seen in 

Fig. m 



V. APPLICATION: THE NATURE OF THE 
INTRINSIC SHUNT 

The Resistively and Capacitivcly Shunted Junction 
(RCSJ) model shows how a sinusoidal supercurrent- 
phase relation, a linear quasiparticle current and a 
displacement-current determine the shape of the entire 
IV characteristic of a Josephson junction. The model 
can also be applied to an unshunted junction, but then 
the subgap resistance R sg appears in the expression for 
the Stewart-McCumber parameter 



(I c Rn) C f R 



Rn 



(50) 



where C is the capacitance of the junction and 'I'oC— 
2.07 x 10~ 15 Wb) the flux quantum. Likhare\^ showed 
that the relation between (3c and the presence of hys- 
teresis depends on the model that is used to describe the 
junction (e.g. Non-linear Resistive model, with different 
dependencies for the subgap conductance, and the Tun- 
nel Junction Microscopic model), but roughly speaking, 
it can be said that (3c > 1 corresponds to hysteretic IV 
characteristics. Hysteresis refers here to the existence of 
two brances in the IV curve, one going from the critical 
current I c to the voltage state, and one going back at the 
return current (Ir < I c ) from the voltage state to the 
state at V = 0. 

The capacitance of a double-barrier junction is not 
known a priori. In Ref. |H a set of Nb/Al double- 
barrier Josephson junctions was fabricated in order to 
make SQUIDs. From resonances in the SQUID washer, 
C was determined to be 0.015 pF/^m 2 , corresponding 
to the capacitance of two SIS junctions in series^ It is 
assumed that this value is only weakly depending on the 
transparency of the barrier. The dependence of I c and 
Rn on the junction parameters, such as jcs, follow from 
the modeling of the stationary properties in Ref. I2H 
The subgap conductance as function of the suppression 
parameter is determined in section Till Bl 

First, the regime of junctions with j c g 3> 1 will be 
discussed. For this purpose, low critical current density 
Nb/Al double-barrier junctions were fabricated accord- 
ing to the process of Ref. Figure EI] shows a typical 
measured IV characteristic, together with an TV curve 
from the nonequilibrium model of section IIVBI where 
inelastic scattering is taken into account. At 4.2 K, the 
experimental and theoretical curves are very much alike, 
taking into account the fact that only one free parameter 
was used to fit, namely Tri n . From 



Ft- — 

- 1 - ' i.n. — 



7eff h 



(51) 



and the fitted TT in =0.1 and 7 ff = 2 x 10 3 (which 
was obtained from fitting the critical current temperature 
dependence), an inelastic scattering time r in = 0.3 ns is 
obtained in the Al interlayers. At 1.6 K, a magnetic field 
was used to suppress the supercurrent in order to resolve 
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the subgap quasiparticle conductance. The deviation of 
the fit from the experiment around 2An6 is due to the 
nonequilibrium enhancement of the gap in the interlayer, 
as described in Ref. which can be included in the 
model by incorporating Aai- However, the good fit well 
below 2Ajvh allows for the extraction of T; n = 0.9 ns at 
1.6 K. 

The values for inelastic scattering correspond to mea- 
surements by Santanam et al~ who found Ti n = 0.2 
to 1.0 ns in 10 nm Al films at 4.2 K, and Van Son et 
(//ii who found Ti n — 0.8 to 0.9 ns in 7 nm Al films at 
T c ai- Our values of Ti„ = 0.9 ns at 1.6 K and 0.3 ns 
at 4.2 K indicate a scaling with T _1 rather than T~ 3 , 
which was found and discussed as well by Santanam et 
alm^ Note, that the values for the inelastic scattering are 
much smaller than 7rfcsT, which means that the station- 
ary properties are not influenced by Ti„. Furthermore, 
from these values it is seen that 7 ff<5~ 1 3> 1 as long as 
7cff < 10 3 , which was used in order to obtain Fig.^jJ 

As a measure of the subgap conductance, the theoret- 
ically expected normalized conductance at eV = Ag can 
be found in Fig. [8] as function of the inelastic scattering 
parameter Tr^ . It can be seen that the subgap resistance 
in the limit of zero inelastic scattering {ji n — > oo) is only 
determined by temperature. The conductance in this 
limit is therefore called the thermal contribution. The 
relation between subgap resistance and 7 e ff is now known 
for a fixed value of Tj„ since T is given by 7rfcsT/7 e ff. 

The dependence of I c Rn on j e g is known from the 
Matsubara modeling of the stationary properties of 
double-barrier junctions. 26 Together with the definition 
of 7eff, it follows that 



R 



N 



e 2 kp Trk B T cS d 



(52) 



where the parameter values can be taken as vp = 1.5 x 



10 6 m/s, 52 d = 6 nm, and T cNb = 9.2 K. Putting these 
theoretical dependencies together with the experimen- 
tally determined parameters into Eq. (|50|l . provides (3c 
as function of the critical current density for junctions 
with 7 off > 1, see Fig.lTTl 

The shunting behavior can physically be explained as 
follows. A direct transfer process of quasiparticles from 
one electrode to the other is prohibited when the quasi- 
particle energy falls within the gap of the other elec- 
trode. However, by scattering inelastically in the inter- 
layer, the quasiparticles are redistributed over energy, al- 
lowing some quasiparticles to enter the other electrode, 
which results in an enhanced conductance, as illustrated 
in Fig. [7| The amount of quasiparticles that get scattered 
inelastically increases for decreasing barrier transparen- 
cies, since the effective lifetime of a quasiparticle in the in- 
terlayer is then increased. For strong inelastic scattering, 
the double-barrier junction can be regarded as a series 
connection of an SIN and NIS junction, where the energy 
distibution function in the interlayer is the equilibrium 
Fermi function /q = tanh(£y2fcsT). Here it should be 
note that the assumption is made that the inelastic scat- 
tering is dominated by electron-phonon interactions, and 
that there is coupling between the interlayer and a heat 
bath. It is known^ that in the contrary non-adiabatic 
limit of MAR and strong electron-electron interactions, 
the energy distribution in the interlayer is given by the 
Fermi function at temperature k B T — A + eV. 

In order to understand the intrinsic shunt of all Al- 
based double-barrier junctions, the regime of high-J c 
junctions (typically larger than 100 A/cm 2 ) should be 
considered as well. The second contribution to the 
subgap conductance is due to the Andreev reflection 
processes at the two superconductor-normal metal in- 
terfaces, which was formally introduced by the term 
Re(i ;l )Re(Fs). The Andreev channels open at high trans- 
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1 10 100 1000 10000 



J c [A/cl/] 

FIG. 12: Theoretical model for /3c as function of J c (and 
as function of 7 e ff in the inset), based on the contribution 
of Andreev channels to the subgap conductance in high-J c 
junctions at T = 4.2 K, in comparison with the hysteresis of 
SIS junctions. 



1 n 


"I ' 


i i mui — i r 


MINI 1 1 1 

. o 


1 ■ ■ ■ ■ | 1 rTTTTTTj 1 /l.'l 1 










o •' !■' ■ 


0.9 






°\ 


0.1 ns / •' ;j p 
/ •*'' / 


0.8 








A '■/■■'' 11 ' 




















\¥ •: □ - 


0.7 








° •'' ■'' / 








0.3 ns'\ 


•''■'/ 


0.6 








,' / 

/ / 




■ T- 

in 


= 1.0 ns 






0.5 








□ ■' / 












0.4 


id i_i 


* 1 l_ 


J-LLLllJ ■ 1 1 


mill * 1111 


1 


10 


100 


1000 10000 



J c [A/cV] 

FIG. 13: Theoretical model (dotted line) for the ratio of 
return and critical current, based on the sum of the shunt- 
ing contributions at 4.2 K from both inelastic scattering and 
Andreev channels, for several inelastic scattering times. Ex- 
perimental data are shown from this paper (■), Ref. [(H] (O)i 
Ref. E3 (□) and Ref. (A). 



parency of the interface barriers. In first order, this con- 
tribution is independent of temperature, but it depends 
on the suppression parameter, which is shown in the in- 
set of Fig. 0] for a fixed temperature. For the practical 
range of parameters, this means that the contribution is 
inversely proportional to 7 c ff. FigureEl s h° ws the result- 
ing hysteresis as function of critical current density. Fig- 
ure 1121 predicts that non-hysteretic double-barrier junc- 
tions can be obtained with critical current densities of 
the order of 10 kA/cm 2 and higher. In order to make 
the comparison with SIS junctions, a similar curve has 
been calculated based on Eq. l(50"|) and plotted in the 
same Fig. ED In this calculation it was assumed that 
C = 3%F/cm 2 , I C R N = 2.0 mV and R sg = 2R N . A big- 
ger subgap resistance will shift the SIS curve even more 
to the right. 

Summing up all contributions to the subgap conduc- 
tance provides the theoretical curve in Fig. ll3l for several 
values of Ti n and d — 6 nm, where Zappe's equation 60 
was used to calculate the ratio of return and critical cur- 
rent Ir/I c from (3c- The summation is performed in a 
straightforward way, since the contributions due to in- 
elastic scattering and the opening of Andreev channels 
do not overlap, i.e. they occur in separate regimes of the 
suppression parameter. 

An increase in Ti n is seen to increase the hysteresis 
and shift the maximum hysteresis to lower values of J c . 
A thicker interlayer will both decrease J c , as well as shift 
the curve upward since Tri n is larger in this case. A 
decrease in temperature rapidly enhances the hysteresis, 
since both the thermal contribution to the subgap con- 
ductance decreases as well as the contribution of inelastic 
scattering, since r; n increases with temperature. This ex- 



plains the strong influence of temperature on hysteresis 
as observed in experiments, which is stronger than could 
be expected from an increase in I c alone. 

Observed experimental Ir/I c values are shown as well 
in Fig. E] and it can be concluded that the experiments 
are now qualitatively and quantitatively very well ex- 
plained by the model in the sense that both the non- 
monotonic hysteresis dependence on critical current den- 
sity as well as the actual hysteresis values are obtained. 



VI. CONCLUSION 

Time-dependent and nonequilibrium transport proper- 
ties of SINIS junctions have been studied by means of a 
microscopic Green's function approach. 

The kinetic equations for the longitudinal and trans- 
verse energy distribution functions are derived from the 
Keldysh-Usadel equation. The appropriate boundary 
conditions are derived by starting from the Kupriyanov- 
Lukichev boundary conditions and by applying Gauge 
transformations in the electrodes. The resulting set of 
equations has a recurrent nature, in terms of coupling 
of Green's functions to higher harmonics as well as to 
functions with shifted energy arguments. This lays out 
a theoretical framework to microscopically study time- 
dependent problems in superconducting - normal metal 
devices. 

We apply this formalism in order to develop a the- 
ory of the subgap conductance of SINIS junctions. This 
conductance is a very favorable feature for applications 
but so far not understood on a microscopic level. In the 
adiabatic limit of a small voltage and a large suppres- 
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sion parameter, the time-dependencies simplify and the 
equations are solved to determine the dissipative current 
in double-barrier Josephson junctions. Known limiting 
cases, such as the SININ' junction, are reproduced. Ex- 
cess and deficit current are determined as function of the 
suppression parameter and the asymmetry between the 
barriers. Excess current as high as eI ex R,N — 1.05Ag can 
exist in double-barrier junctions in the symmetric case for 
7eff <C 1, and maximum deficit current is reached in the 
symmetric case for %g ^> 1. The subgap conductance 
enhancement by decreasing j e g is caused by the opening 
of Andreev channels. 

It is found that the time-dependent nonequilibrium 
contribution to the energy distribution function gives rise 
to a nonzero averaged supercurrent in the presence of a 
voltage bias. This effect should be observable in double- 
barrier junction experiments. 

In contrast to most studied mesoscopic systems, inelas- 
tic scattering in the interlayer of double-barrier junctions 
can have a strong influence on the electronic transport 
even in very short devices. A microscopic derivation of 
the dependence of the transport properties on the inelas- 
tic scattering parameter is given. IV characteristics show 
an enhanced subgap conductance for increased inelastic 
scattering rates. 

The actual value for the inelastic scattering time in 
the interlayer of experimentally realized devices was ob- 
tained by fitting the microscopic model. The inelastic 
scattering values explain, together with the opening of 
Andreev channels, the nature of the intrinsic shunt in 
double-barrier junctions. 
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APPENDIX A: DERIVATION OF THE 
BOUNDARY CONDITIONS 

From the definition of the Green's functions in Keldysh 
x Nambu space, Eq. l(T)l. a boundary condition can be 
written for each of the matrix elements of Eq. 124|l 

lB 0n^G R = G R G R ~ G R G R , 
ax 

ax 

7B £ + 

\ ax ax 
= G fl Gf + G K G A - G R G K - Gf G A . (Al) 



With definition (|18fl . the left-hand side of the latter of 
these three boundary conditions becomes 



G RjL G x)f + G R G R -^f 
ax I ax 



( f G A + G R f^-G A + fG A -^-G J 
ax V / ax ax 



(A2) 



With the aid of the first two Eqs. in (|A1|) this can be 
rewritten into 



Jb£, 



G n G Rdj_ G Rd 
ax ax 



(fG A ) 



-G R f ——G A 

ax 



+ (d R G R - Gf G fl ) / - / (d A G A - G A G A ^j .(A3) 

By making use of the normalization condition 
G R(A) G R(A) = x and the definition for / an d f u 

Eq. (|19f) . this can be futher rewritten as 



1b£ 



^ (ft + r 3 fr) - G R ^ (f L + f 3 f T ) G A 



\G R G R -G R G R ) (h+r 3 f T ) 
-(/l+t 3 / t ) (g a G a -G a G a 
which is equal to 



(A4) 



1b£ 



(l - G R G A ) ±f L + (r 3 - G R f 3 G A ) ±f 7 



GR./^iR. r-ili.r~ili. r~ij\r~ij\ . r^Ar^A \ f 
Lt-l ' ' | \ ' y i Lr-^ t Lr-^ Lr I J L 



'<Rr<R 



AAA 



^iAAA 



G R G R - Gf G«) f 3 - f 3 (G A G A ~ G A G A 



It- 
(A5) 



The right-hand side of the last boundary condition in Eq. 
(IA1I) can be rewritten with the aid of Eq. i|18|) into 



G R G R A - Gf f x G A - G R hG A + fcG A G A 



fG A G A - Gf G n f + G"/Gf + Gf /G 



R ifiA 



^R triA 



(A6) 



TLl 



Tti 



With definition Eq. (|19|l for / and /i this becomes 
G R (Gf - Gf ) - (Gf - Gf) G A ] I 

G R (Gf f 3 - r 3 G A ) (Gf f 3 - f 3 Gf) & 
(G R - G A ) Gf - Gf (G R - G A ) 

(A7) 

Equating left-hand and right-hand side of the last bound- 
ary condition in Eq. (jAlfl . i.e. Eqs. (|A5f) and (|A7|) 
respectively, finally gives the form of the boundary con- 
dition as presented in Eq. I|25|) . 
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APPENDIX B: TIME CONVOLUTIONS IN 
ENERGY SPACE 



With an energy-shift E = E + uujq/2 this becomes 



The expression for a convolution of two functions, 

oo 

aob(t,t')= / diia(i,ti)6(ti,i'), (Bl) 



can be transformed by changing variables 
aob(t,t') = 



,// L „ ( / . / 1; LJl) b ( tl - t', t -^- } . (B2) 



Subsequently, a Fourier transform to energy-frequency 
space can be made 



aob(t, t') = J dtxdujdEdjj'dE'a (E, w) e lE{t - t ^l h 

— OO 

e iu(t+ti)/2K b ^ ^ e iE'(t 1 -t')/h e i "'(* 2 1 h +t '') 

OO 

= E/ dhdEdE'an (E) e mt-ti)/H e inu, (t+ tl )/2H 
b i (E' s j e iE '{ t i- t ')/ h e in ' ul °{ t i-+ t ')/ 2h 

OO 

= E / dEdE ' a n (E) e iEt/h e inw t/m bni ^ 
"'"'-00 

e -iE't'/H^n'w t'/2H 5 ^_ E + nujQ / 2 + E l + n ' Wo / 2 ) 

OO 

= E / dE ' a - ( E ' + bn ' 

"'"'-00 

e ~iE\t~t')/h e i^u {t-t')/h e i^^u {t+t') 



OO 

oo6(t,0 = E/ dEa n (E> + ^y n ,(E-^ 



xe 



) 

(B4) 



(B3) 



The triple products in the boundary conditions (3.34) 
and (3.35) can be worked out in the same manner. The 
sine and cosine terms cause shifts in the arguments. An 
example of the result of a triple convolution is given in 
Eq. (ED. 
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